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Abstract 

In  this  paper  we  focus  on  the  long-term  behavior  of  generalized  polynomial  chaos 
(gPC)  and  multi-element  generalized  polynomial  chaos  (ME-gPC)  for  partial  differ¬ 
ential  equations  with  stochastic  coefficients.  First,  we  consider  the  one-dimensional 
advection  equation  with  a  uniform  random  transport  velocity  and  derive  error  esti¬ 
mates  for  gPC  and  ME-gPC  discretizations.  Subsequently,  we  extend  these  results 
to  other  random  distributions  and  high-dinrensional  random  inputs  with  numerical 
verification  using  the  algebraic  convergence  rate  of  ME-gPC.  Finally,  we  apply  our 
results  to  noisy  flow  past  a  stationary  circular  cylinder.  Simulation  results  demon¬ 
strate  that  ME-gPC  is  effective  in  improving  the  accuracy  of  gPC  for  a  long-term 
integration  whereas  high-order  gPC  cannot  capture  the  correct  asymptotic  behav¬ 
ior. 
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1  Introduction 


Polynomial  chaos  (PC)  has  been  used  extensively  in  the  last  decade  to  model 
uncertainty  in  physical  applications  [1-6].  It  is  based  on  the  original  ideas  of 
homogenous  chaos  (Wiener- chaos)  first  formulated  by  Wiener  as  the  span  of 
Hcrmite  polynomial  functionals  of  a  Gaussian  process  [7].  Ghanem  &  Spanos 
were  the  first  to  combine  Wiener-chaos  with  a  finite  element  method  to  model 
uncertainty  addressing  solid  mechanics  applications  [1,8,9].  A  more  general 
framework,  termed  generalized  polynomial  chaos  (gPC),  was  proposed  in  [10] 
by  Xiu  &  Karniadakis,  following  the  framework  of  Ghanem  &  Spanos,  based 
on  the  correspondence  between  the  PDFs  of  certain  random  variables  and  the 
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weight  functions  of  orthogonal  polynomials  of  the  Askey  scheme.  The  family 
of  gPC  includes  Hcrmite-chaos  as  a  subset  and  provides  optimal  bases  for 
stochastic  processes  represented  by  random  variables  of  commonly  used  dis¬ 
tributions,  such  as  uniform  distribution,  Beta  distribution,  etc.  Polynomial 
chaos  was  combined  with  wavelets  in  [11,12]  to  deal  with  discontinuities  for 
uniform  random  inputs  for  which  standard  PC  or  gPC  fails  to  converge.  To 
solve  differential  equations  with  stochastic  inputs  following  the  procedure  es¬ 
tablished  by  Ghanem  &  Spanos,  the  random  solution  is  expanded  spectrally 
by  polynomial  chaos  and  a  Galerkin  projection  scheme  is  subsequently  used 
to  transform  the  original  stochastic  problem  into  a  deterministic  one  with  a 
large  dimensional  parameter  [1,10,6]. 

On  the  other  hand,  Deb,  Babuska  &  Oden  [13]  have  proposed  to  employ  finite 
elements  in  the  random  space  to  approximate  the  stochastic  dependence  of 
the  solution.  This  approach  also  reduces  a  stochastic  differential  equation  to 
a  high  dimensional  deterministic  one.  This  method  was  later  studied  theo¬ 
retically  within  the  framework  of  deterministic  finite  element  method  in  [14]. 
Since  a  finite  element  method  is  generally  used  to  solve  the  obtained  deter¬ 
ministic  PDE  system,  the  above  methods  are  called  stochastic  Galerkin  finite 
element  method  in  [14]  while  the  scheme  in  [1,10]  is  classified  as  p  x  h  version 
and  the  scheme  in  [13]  as  k  x  h  version.  Here  p  denotes  the  polynomial  or¬ 
der  of  polynomial  chaos,  k  the  element  size  in  the  random  space,  and  h  the 
element  size  in  the  physical  space.  Both  schemes  use  finite  elements  in  the 
physical  space.  The  px  h  version  relics  on  the  global  representation  in  the  en¬ 
tire  random  space  by  polynomial  chaos  while  the  k  x  h  version  is  based  on  the 
discretization  of  the  random  space  using  the  same  basis  as  the  deterministic 
finite  element  method  to  approximate  the  random  held  locally.  Both  concepts 
and  the  terminology  introduced  here  have  similarities  with  the  spectral / hp 
element  method  for  deterministic  problems  [15,16]. 

Although  gPC  works  effectively  for  many  problems,  e.g.,  elliptic  and  parabolic 
PDEs  with  stochastic  coefficients  [14,6],  it  cannot  deal  with  some  other  dif¬ 
ferential  equations,  e.g.,  the  Kraichnan-Orszag’s  three-mode  ODE  system  for 
modeling  turbulence  [17]  or  the  Navier-Stokes  equations  for  unsteady  noisy 
hows  such  as  how  past  a  stationary  cylinder  [18].  For  these  problems  gPC 
fails  to  converge  after  a  short  time,  and  increasing  the  polynomial  order  helps 
little  for  the  convergence.  There  exist  at  least  two  reasons  for  the  divergence 
of  gPC: 

(1)  Singularity  in  the  random  space,  and 

(2)  Long-term  integration. 

The  former  was  studied  in  [11,12]  with  Wiener-Haar  expansions  and  in  [19,20] 
with  ME-gPC.  In  this  work,  we  focus  on  the  second  one:  long-term  integra¬ 
tion.  In  particular,  we  are  interested  in  the  cases  which  are  related  to  random 
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frequencies.  We  use  a  simple  one-dimensional  advection  equation  with  a  uni¬ 
form  random  transport  velocity  as  a  model  problem.  We  first  derive  the  error 
estimates  of  gPC  and  ME-gPC  for  the  Legendre-chaos  expansion.  Based  on 
these  error  estimates  we  obtain  a  relation  between  gPC  and  ME-gPC,  which 
we  verify  by  numerical  computations.  We  subsequently  generalize  such  a  re¬ 
lation  to  other  random  distributions  and  more  general  random  inputs.  Lastly, 
we  consider  a  physical  problem:  noisy  flow  past  a  stationary  circular  cylinder, 
where  we  use  the  obtained  results  to  analyze  the  convergence  of  gPC  and 
ME-gPC. 

This  paper  is  structured  as  follows.  In  section  2  we  present  an  overview  of 
gPC  and  ME-gPC.  In  section  3  we  study  theoretically  the  one-dimensional 
stochastic  advection  equation.  In  section  4  we  present  results  for  stochastic 
simulations  of  noisy  flow  past  a  stationary  circular  cylinder.  We  conclude  with 
a  short  discussion  in  section  5. 


2  Overview  of  gPC  and  ME-gPC 


Let  (f2,  T ,  P)  be  a  complete  probability  space,  where  is  the  sample  space,  T 
is  the  ex-algebra  of  subsets  of  and  P  is  a  probability  measure.  An  Md-Valued 
random  variable  is  defined  as 

Y  =  (WM, . .  .,Yd(u))  :  (n,p)  »  ( Rd,Bd ),  (1) 

where  d  €  N  and  Bd  is  the  u-algebra  of  Borel  subsets  of  Md.  If  uncertainty 
is  included  in  the  PDEs,  the  solutions  can  be  referred  to  as  a  random  field 
u(x,  t‘,u),  where  x  denotes  the  physical  space  and  t  the  time.  For  any  fixed  x 
and  t,  u(x,  t ;  u>)  is  a  Mdp-valued  random  variable,  where  dp  is  the  dimension  of 
physical  domain  (dp  <  3).  We  generally  assume  that  u(x,  f;  u)  is  a  second-order 
random  field  denoted  as  u(x,  f;cu)  G  P,  P),  where 

[  u2(x,  t]  u))dP(uj)  <  oo.  (2) 

Ju 

We  expect  that  all  viscous  flow  satisfy  this  constraint. 


2.1  Generalized  polynomial  chaos  (gPC) 


Generalized  polynomial  chaos  is  a  spectral  polynomial  expansion  for  a  second- 
order  random  field  u(x,  t]u),  i.e., 

OO 

u(x,  t;uj)  =  J2  a*(x,  (a;)),  (3) 

i=0 
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where  {<£>j(W)}  denote  the  basis  of  gPC  in  terms  of  Y ;  this  spectral  expansion 
converges  in  the  L2  sense.  We  usually  select  a  weighted  orthogonal  system 
{<£>j(W)}  in  L'2 (hi,  JF,  P )  satisfying  the  orthogonality  relation 

(4) 

where  St]  is  the  Kronecker  delta,  and  (•,  •)  denotes  the  ensemble  average  with 
respect  to  the  probability  measure  P.  The  index  in  equation  (3)  and  d  G  N 
are,  in  general,  infinite.  In  practice,  both  limits  will  be  truncated  at  a  certain 
level. 

For  a  certain  Md-valued  random  variable  Y ,  the  gPC  basis  {<£>*}  can  be  chosen 
in  such  a  way  that  its  weight  function  has  the  same  form  as  the  probability 
density  function  (PDF)  of  Y.  The  corresponding  type  of  classical  orthogonal 
polynomials  {<£>,}  and  their  associated  random  variable  Y  are  listed  in  table  1 
[10].  For  arbitrary  probability  measures,  the  orthogonality  must  be  maintained 
numerically,  as  we  explain  in  the  next  subsection. 


Table  1 

Correspondence  of  the  type  of  Wiener- Askey  polynomial  chaos  and  their  underlying 
random  variables  (TV  >  0  is  a  finite  integer). 


Random  variables  Y 

Wiener- Askey  chaos  {^(l^)} 

Support 

Continuous 

Gaussian 

Hermite-chaos 

(—00, 00) 

gamma 

Laguerre-chaos 

[0, 00) 

beta 

Jacobi-chaos 

[a,b\ 

uniform 

Legendre-chaos 

[a,b\ 

Discrete 

Poisson 

Charlier-chaos 

{0,1,2,...} 

binomial 

Krawtchouk-chaos 

{0,1,..., TV} 

negative  binomial 

Meixner-chaos 

{0,1,2,...} 

hypergeometric 

Hahn-chaos 

{0,1,..., TV} 

2.2  Multi- element  generalized  polynomial  chaos  (ME-gPC) 


We  assume  that  Y  is  defined  on  B  =  xf=1[a.j,  bj],  where  a*  and  bi  are  finite  or 
infinite  in  M  and  the  components  of  Y  are  independent  identically -distributed 
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(i.i.d.)  random  variables.  We  define  a  decomposition  D  of  B  as 


D 


Bk  —  X  [Ofc,2,&fc,2)  x  •••  X  [dk,d,bk,d], 


Bk]  fl  Bk2  =  0,  if  hi  7^  k2, 


(5) 


where  k,ki,k2  =  1,2,  •••  ,  IV.  Based  on  the  decomposition  D,  we  define  the 
following  indicator  random  variables 


lB, 


1  if  YeBk, 

0  otherwise. 


(6) 


Thus,  SI  =  U iLjBl  (1)  is  a  decomposition  of  the  sample  space  O,  where 


ib](X)  n  ib](X)  =  0,  for 


(7) 


Subsequently,  we  define  a  new  Revalued  random  variables  C,k  :  /s  "(1)  - 
Bk  on  the  probability  space  (/s*( 1),.F  fl  (1) ,  -7)(*  1 1  Bk  1))  subject  to  a 
conditional  PDF 


/fc(y|7Bfc  —  i)  - 


f(y ) 


(8) 


Pr(/Bfe  =  1)  f 

where  f(y )  denotes  the  PDF  of  Y  and  Pr (IBk  =  1)  >  0.  In  practice,  we  usually 
map  Ck  t°  a  new  random  variable  Y k  defined  on  [—1,  l]d  to  avoid  numerical 
overflow  in  computer  [20],  using  the  following  linear  transform 


Ci  =  g(Yi) :  Ci 


bk,i  dkyi-y  bk,i  ~f" 


(9) 


where  i  =  1,  2,  •  •  •  ,  d  and  k  —  1,  2,  •  •  •  ,  N.  To  this  end,  we  present  a  decom¬ 
position  of  random  space,  which  is  very  similar  with  the  decomposition  of 
physical  space  using  separable  elements. 


Based  on  the  random  variables  {W^},  a  scheme,  called  multi-element  general¬ 
ized  polynomial  chaos  (ME-gPC),  was  proposed  in  [19,20].  Based  on  ME-gPC, 
u(x,  t]  u>)  can  be  expressed  as  [19] 

N  M 

u(x,  t;  u>)  =  ^2^2dkti(x,t)$kji(Yk(Y))IBk,  (10) 

k= 1  i=0 

where  {<!>*.;}  is  the  local  chaos  basis  in  element  k  and  M  is  the  number  of 
chaos  modes.  The  key  idea  of  ME-gPC  is  to  implement  gPC  element-by- 
element  when  the  global  spectral  expansion  is  not  efficient  to  capture  the 
random  behavior.  Thus,  the  basic  procedure  of  ME-gPC  is  quite  similar  with 
the  deterministic  spectral  element  method.  In  the  decomposition  of  physical 
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space  using  continuous  Galerkin  projections,  we  need  to  treat  carefully  the 
connectivity  (C°  continuity)  between  two  adjacent  elements;  however,  in  the 
decomposition  of  random  space  the  following  C'°-type  continuity 

ufll(r)  =  ufl2(r),  YeB,  db2,  (li) 

where  Bj  is  the  closure  of  element  BtJ  is  not  required  since  the  Lebesgue  mea¬ 
sure  of  the  interface  between  two  random  elements  is  zero  and  most  statistics 
we  are  interested  in  are  defined  as  a  Legesgue  integration. 

In  the  decomposition  of  random  space,  the  PDF  of  Y  is  decomposed  simulta¬ 
neously,  which  implies  that  the  original  gPC  basis  will,  in  general,  lose  local 
orthogonality  in  random  elements.  The  only  exception  is  the  Legendre-chaos 
for  the  uniform  distribution  [19].  For  other  distributions,  orthogonal  polyno¬ 
mials  with  respect  to  the  PDF  of  local  random  variable  Y *.  can  be  constructed 
numerically.  Given  an  arbitrary  PDF,  the  Stieltjes  procedure  and  the  Lanczos 
algorithm  [21]  can  be  used  to  construct  the  following  orthogonal  system 


TTi+i(t)  =  (t  -  ai)7Ti(t)  -  PiTTi-^t),  i  =  0, 1,  •  •  •  , 

n0(t)  =  l,  vr_1(t)  =  0,  (12) 

where  {^(f)}  is  a  set  of  (monic)  orthogonal  polynomials, 

7 Ti(t)  =  tf  +  lower-degree  terms,  i  —  0, 1,  •  •  •  (13) 

and  the  coefficients  cq  and  j3i  are  uniquely  determined  by  a  positive  (prob¬ 
ability)  measure.  The  orthogonal  system  {7q (Y")}  will  serve  as  a  local  gPC 
basis. 

In  ME-gPC,  relative  low  polynomial  orders  (5  to  8)  are  preferred  locally;  thus, 
the  numerical  re-construction  can  be  implemented  efficiently  and  accurately 
[20].  Numerical  experiments  show  that  the  cost  of  maintaining  local  orthogo¬ 
nality  is  negligible  compared  to  the  cost  of  a  standard  Galerkin  gPC  solver. 


3  Long-term  integration  of  gPC  and  ME-gPC 


In  this  section  we  study  the  long-term  behavior  of  gPC  and  ME-gPC  analyt¬ 
ically  using  the  following  one-dimensional  stochastic  advection  equation 


du 

c)t 


+  no 


du 

dx 


0 


subject  to  the  initial  condition 


u(x)  =  u0(x-,^), 


(14) 


(15) 
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where  £  is  a  one-dimensional  uniform  random  variable  defined  on  [—1,1]  and 
P(£)  G  I/2 (0,  P,  P).  In  particular,  we  assume  that 

V(£)=u  +  o-£,  u0{x]£)  =  sinn7r(l  +x),  a:  G  [-1, 1],  (16) 

where  a  is  a  constant,  v  is  the  mean  of  transport  velocity  and  n  e  N.  It  is 
easy  to  obtain  the  exact  solution  for  this  case  as 

u(x,  t]  £)  =  sinn7r(l  +  x  —  (h  +  cr£)t),  (17) 

which  shows  that  the  frequency  of  this  stochastic  process  is  random. 


3. 1  Error  estimates  for  gPC 


Let  (Pj(£)}  denote  the  orthogonal  basis  of  Legendre-chaos  and  Vm  denote  the 
projection  operator  as 


M 

VMu{x,  t-  £)  =  ui(x >  (18) 

i=0 

where 

=  7p^|u /1M(x>i;^)-p<(^)^-  (19) 

We  study  the  convergence  of  VMu(x,t]£)  since  it  can  demonstrate  the  main 
properties  of  numerical  convergence. 


Theorem  1  Let  6m  denote  the  error  of  the  second-order  moment  of  Vmu. 
Given  time  t  and  polynomial  order  M,  €m  can  be  bounded  as 


eM  <  C(M) 


2M+2 
HM 

1  ~<lV 

where  C(M )  is  a  constant  depending  on  M  and 

aniret 


(20) 


Qm  = 


2  M  +  2 

Here  e  is  the  base  of  natural  logarithm. 


<  1. 


(21) 


PROOF.  According  to  the  following  formula  in  [22] 


sinc7r(z  +  a)  =  — p=  ^(2 i  +  1)  A+i/2(c7t)  sin(c7ra  +  -m)Pi(z),  (22) 

v 2c  i=0  2 
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we  obtain  the  polynomial  chaos  expansion  of  the  exact  solution  (17)  as 


OO 

“=-yssg(2l+ 

where  Ji+1/2  are  Bessel  functions  of  the  first  kind.  Using  the  orthogonality  of 
Legendre  polynomials,  we  obtain 

1  00  1 

(u2(x,  t;  0)  =  +  l)Ji+ i/iianirt)  sin2(mr(ut  -  x  -  1)  +  -in),  (24) 

where  (Pi(C)Pj(0)  =  / (27  +  1)  is  employed.  Then  6m  can  be  expressed  as 


1)  Ji+i/2(cmnt)  sm(nn(vt  -  x  -  1)  +  -in)Pi(£),  (23) 


6m  =  {u  )  —  (( Vmu )  ) 

1  00  1 

y,  (2i  +  1)  Jf+1/2(annt)  sm2(nn(vt  —  x  —  1)  4 — in).  (25) 


2 nat  i=M+i 

It  is  known  (see  [23])  that 


2anntJi+1/2^an7Tt>>  =  2^i!  L  cos((crn7ri)  cos(#))  sin2*+1  6>d6>.  (26) 

By  substituting  equation  (26)  into  equation  (25),  cm  can  be  approximated  as 

^  (2i  +  l)(annt)2\  .  2  ,  1 

eM=  y  22,+2(?;,)2 - Asm  (nn(vt^.x-l)  +  -m), 


i=M+ 1 


Ai=i  cos(annt  cos(#))  sin2l+1  OdO  I  . 


where 


Using  Stirling’s  formula  [23]  for  the  factorial  ?'!,  we  obtain  that 
^  (2i  +  l)(crnnte)21  .  2  ,  1 

i.tr+1  8^(2 if  '  '  2  '• 

where  e  is  the  base  of  natural  logarithm.  For  a  fixed  time  t ,  the  error  cm  can 
be  bounded  as 

“  (2M  +  3)(qn7ref)2*  (2M  +  3)g^/+2 

-  1  /  J  0  /  7,  ,T  ,  -|  \  71  ,T  I  C\\r)^  1 


i=M+ 1 


8tt(M  +  1)(2M  +  2)2i  ~L%n(M  +  1)(1  -  q2M)  ’ 


where  Ci  is  a  constant  and  qM  =  crnnet/(2M  +  2).  Here  the  condition  qm  <  1 
is  assumed  for  the  convergence  of  summation.  We  subsequently  check  the 
constant  C\ .  Since  sin  0  >  0  in  0  e  [0,  n\,  we  obtain  that 

A1/2  <  r  sin2*+1  OdO. 

Jo 
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Let  Bi  =  Jq  sin2?+1  9d9.  Using  sin2  9  +  cos 2  9  =  1,  the  following  relationship 
can  be  obtained 

Bi  =  sin2*-1  9d9  -  r  sin2*"1  6  cos2  9d9  =  Bt_ i  -  T  sin2*"1  9  cos2  0dd. 

Jo  Jo  Jo 

Since  the  second  term  on  the  right-hand  side  is  positive,  we  know  that  the 
sequence  {Bi}  is  decreasing.  Thus,  we  can  bound  Ai  as 

A<B2<B2M+1,  i  >  (M  +  1). 

Let  Ci  =  B\,1+x  and  C(M)  =  Ci(2M  +  3) / 87 r(M  +  1),  then  the  conclusion 
follows  immediately. 


In  Theorem  1,  qm  <  1  is  assumed  for  the  convergence  of  summation  in  equa¬ 
tion  (3.1).  For  a  general  case  we  have  the  following  corollary: 

Corollary  2  Given  time  t  and  polynomial  order  M ,  6m  can  be  bounded  as 

1  m  „  q2fI+ 2 

evr  <  —  Y,  ( 2i  +  +  C(M)  ^  2  ,  (27) 

Z(jr  j=M+l  1 


where  is  a  function  of  M  defined  as  in  equation  (21)  and  q ^  <  1. 


3.2  Error  estimates  for  ME-gPC 


Let  Vm  denote  the  projection  of  u(x,t',0  onto  the  basis  of  ME-gPC. 


Theorem  3  Given  a  decomposition  of  random  space  of  f  with  element  length 
Lk  =  bk  —  ak,  k  =  1,  2,  •  •  •  ,  N,  the  error  cm  of  the  second-order  moment  of 
Vm'u  can  be  bounded  as 


N 


eM  <  C(M)  Y 


2M+2 

9k, M 


=1 1  -  %, 


p  t{Ib„  =  1), 


M 


(28) 


where  C(M )  is  a  constant  depending  on  M  and 

crnneLkt 

qkM  =  2(2 M  +  2)  < 


(29) 


PROOF.  According  to  equation  (10),  we  know  that  Vm  can  be  expressed  as 

N 

Vm  =  Y,  'Pk,MlBk ,  (30) 

fc=i 
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where  Vk,M  is  a  local  projection  operator  defined  as 


> \Mu{x,  t-  4)  =  VMU{X,  t ;  — - — 4  H - - — , 

Then,  the  second-order  moment  can  be  expressed  as 


(31) 


N 


({VMu(x,  f;  4))2)  =  (  5Z  t ;  £)IBk  ) 


\k= 1 


AT  ^  2 

=  ;£))  )  Pr(^st  =  i). 


(32) 


A=1 


Thus,  £m  takes  the  following  form 


N 


Uu  —  X]  P^-Tb*.  —  1), 


(33) 


k= 1 


where  a/  is  the  error  of  the  second-order  moment  of  Vk,Mu(x,t]^).  We  now 
check  the  behavior  of  ek,M-  Given  a  random  element  Bk  =  [ak,bk],  the  lo¬ 
cal  problem  in  ME-gPC  is  to  find  the  solution  of  the  following  transformed 
equation 

du  (_  ( bk  —  ak  bk  +  ak\\du 

¥  +  q  +  <7(— >‘  +  — — 

where  4a-  is  the  local  uniform  random  variable  defined  on  [—1,1],  The  exact 
solution  of  equation  (3.2)  is 


u(x,  t ;  4a)  =  sin  m r 


i  i  ( -  i  ( ak  c  .  bk  T  ak  .  .  , 
1  +  x  -  \  v  +  a\  — 4fc  -I - - -  |  |  t 


2  2 

/  / 

Using  a  similar  procedure  as  in  section  (3.1),  we  can  bound  ek>M  as 

C{M)&2 


tk,M  < 


where 


Qk,M  — 


i-  Qk  ,M 

anne(bk  —  ak)t 


<  1. 


2(2  M  +  2) 

Using  equation  (33),  the  conclusion  follows  immediately. 


3.3  Relation  between  6m  and  £m 


From  theorem  1,  we  can  see  that  Qm  increases  linearly  in  terms  of  time  t, 
which  implies  that  gPC  will  lose  /^-convergence  after  a  finite  time.  To  keep  a 
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certain  accuracy,  the  polynomial  order  of  gPC  must  increase  with  time.  Let 


V  C(M)  2M+2 

0  i  „2  "M  i 

1  Qm 


(34) 


where  <5  denotes  a  desired  accuracy.  By  solving  such  an  equation  we  obtain 
that 


'5(1  —  a 2  i\  V(2M+2) 

(2M  +  2). 


t  = 


arnre  y  C(M )  y 
Since  [<5(1  —  q\{) /C(M)]1^2M+2'>  — ■>  1  when  M  — ■>  oo,  we  obtain  that 


2  M  +  2 


arnre 


(35) 


(36) 


which  is  a  linear  relation.  It  is  instructive  to  define  the  increasing  speed  of 
polynomial  order  as 

d  M  arnre 
_  _ 

dt  2 

which  shows  that  to  maintain  an  accuracy  5  the  polynomial  order  must  in¬ 
crease  at  a  speed  arnre/ 2;  we  note  that  n/2  is  the  wave  number  in  the  initial 
condition.  We  can  see  that  the  speed  is  proportional  to  the  wave  number  and 
the  degree  of  perturbation,  which  implies  that  gPC  will  quickly  fail  to  con¬ 
verge  for  a  problem  with  a  large  perturbation  or  wave  number  if  a  random 
frequency  in  time  is  involved. 


Theorem  4  To  maintain  a  certain  accuracy  of  the  second-order  moment  of 
Vmu,  the  polynomial  order  of  gPC  must  increase  with  time  and  the  following 
relation  is  satisfied 


M  «  -  armed  —  1. 
2 


(38) 


We  assume  that  a  uniform  mesh  is  employed  and  p-convergence  is  maintained, 
in  other  words,  qM  <  1  and  qk,M  <  1  with  k  =  1,  2,  •  •  •  ,  N.  Thus,  we  have 

»2M+2 

eM  =  C{M)^r.  (39) 

1  Qm 

The  ratio  of  eM  and  cm,  for  a  fixed  time  t  and  polynomial  order  M,  is 


6m_(1\2M+2  1  ~q2M  /1\2M+2 

eM  \NJ  1  -(jj)2q2M~ 


(40) 


which  is  consistent  with  the  ^-convergence  (e  oc  _/\r-2(M+1))  Qf  ME-gPC  (see 
[13,14,19,20]).  Let  us  consider  that  gPC  and  ME-gPC  of  polynomial  order  M 
reach  accuracy  of  the  same  order.  To  satisfy  this,  we  need  to  have 


Qm  —  Qm, 


(41) 
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which  yields 


tM  —  NtM, 


where  tM  and  tM  denote  time  for  ME-gPC  and  gPC,  respectively. 


(42) 


Theorem  5  Suppose  that  the  error,  e,  of  the  second-order  moment  of  gPC  of 
order  M  is  maintained  in  the  range  t  <tg.  Based  on  a  uniform  mesh  with  N 
random  elements,  ME-gPC  of  order  M  can  maintain  the  accuracy  0(e)  in  the 
range  t  <  Ntg.  In  other  words,  ME-gPC  can  extend  the  valid  integration  time 
of  gPC  linearly  by  a  factor  N. 


If  the  mesh  is  non-uniform,  the  aforementioned  linearity  is  still  valid;  however, 
the  factor  will  be  less  than  N.  We  assume  that  U]f=lBf,.  is  a  decomposition 
for  f  of  uniform  distribution  and  the  length  of  Bk  is  an  increasing  series, 
0  <  l Bi  <  l b2  <  '  "  <  Ibn-  From  the  proof  of  theorem  3,  we  know  that 


N 

=  E  Pr(^Bfc 

k= 1 


N 


1)  <  C(M)  y 

k= 1 


1 


2M+2 


Qk,M 


(43) 


where 


Qk,M 


crmrelBkt 
2(2  M  +  2) 


and  Pr (IBk 


1) 


Ib 


k 


2  ' 


We  define  a  function 


where 


It  is  easy  to  verify  that  Q(z) 
Zk  =  lBk/ 2-  We  can  obtain 


Q{z)  = 


n2M+2 

liz _ 


arnret 


Qz  (2M  +  2) Z' 

is  an  increasing  function  with  respect  to  z.  Let 


N 

Q(z l)  =  Y,  Q(Zl)Zk 


k= 1 


N  N 

<  E  Q(Zk)Zk  C  ^  )  Q(zN)zk  Q(zpj), 

k= 1  k= 1 


where 

N 

YZk  =  1,  0  <  Zi  <  Z2  <  ■  ■  ■  <  ZN. 

k=  1 

To  satisfy  cm  =  £m,  we  need 


Q(z\)  <  Q(  1)  <  Q(zN), 


which  implies  that 

2  „  2 
- - tM  <  tM  <  - tM- 

“Bn  “  B\ 


(44) 
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3-4  Other  distributions  and  high- dimensional  random  inputs 


In  ME-gPC,  the  PDF  of  £  will  be  decomposed  simultaneously  with  the  random 
space;  thus,  the  local  orthogonality  has  to  be  maintained  numerically.  The  only 
exception  is  the  Legendre-chaos  [19]  due  to  the  nice  properties  of  uniform 
distribution.  It  is,  in  general,  difficult  to  analyze  theoretically  the  convergence 
for  the  numerical  basis  of  ME-gPC.  In  this  work,  we  compare  the  performance 
of  gPC  and  ME-gPC  numerically  for  other  distributions. 

The  fc-type  convergence  was  shown  theoretically  in  [14,13]  to  be 

||E[«]  -  E[um}\\l2{d)  <  Ck2(M+1\  ||E[u2]  -  E[u2m]\\L2{D)  <  Ck2^  (45) 

using  an  stochastic  elliptic  model  problem,  where  C  is  a  constant  depending  on 
M,  and  k  denotes  the  maximum  size  of  random  elements.  We  note  here  that  D 
indicates  the  physical  space.  It  was  shown  in  [19,20]  that  the  index  of  algebraic 
convergence  of  ME-gPC  for  the  mean  and  variance  goes  asymptotically  to 
2 (M  +  1)  for  a  uniform  mesh,  which  is  consistent  with  equation  (45).  Note  that 
the  error  bound  (45)  is  independent  of  probability  measures.  Such  observations 
imply  that  for  any  probability  measure  the  following  relation  (see  equation 
(40)) 

—  ~C(M)(4)2M+2  (46) 

6m  N 

holds  for  the  solution  of  equation  (17),  where  the  constant  C  depends  on  the 
polynomial  order  M .  Here  we  include  the  time  t  together  with  the  number 
N  of  random  elements  because  in  the  k-th  random  element  of  ME-gPC  the 
solution  takes  the  form 

u(x,t',Ck)  =  sinn7r(l  +  x  -  (vk  +  (47) 

where  t  and  N  can  be  treated  together.  From  equation  (46)  we  can  see  that 
if  we  re-scale  the  time  of  gPC  by  the  number  N,  6m /^m  will  be  a  constant 
depending  on  the  polynomial  order.  However,  such  a  relation  will  be  reached 
asymptotically  because  of  the  non-uniform  random  distribution  [20]. 

It  is  easy  to  generalize  the  obtained  results  to  high- dimensional  random  inputs. 
Since  the  high-dimensional  basis  of  gPC  is  constructed  by  tensor  products  of 
one-dimensional  basis,  the  error  of  chaos  expansion  should  be  dominated  by 
the  summation  of  errors  of  one-dimensional  truncation.  Thus,  the  results  for 
the  one- dimensional  case  should  be  still  valid  for  a  high- dimensional  case.  For 
example,  if  we  have  d-dimensional  uniform  random  inputs,  ME-gPC  with  Nd 
uniform  elements  should  extend  the  valid  integration  time  of  gPC  with  the 
same  polynomial  order  by  a  factor  N.  However,  in  practice,  the  degree  of 
perturbation  in  each  random  dimension  is  generally  different,  and  only  the 
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random  dimensions  with  large  perturbations  are  needed  to  be  refined.  Such 
cases  are  more  difficult  to  analyze  and  beyond  the  scope  of  this  paper. 


3. 5  Initial  conditions 


For  the  initial  condition  uq(x',£),  we  intentionally  employed  functions  such 
as  cosine  and  sine  waves,  which  introduce  “random  periodicity”  in  time  for 
a  given  random  transport  velocity.  Such  solutions  are  often  encountered  in 
practice,  e.g.,  random  oscillators  and  simulations  of  unsteady  turbulent  or 
noisy  flows.  If  the  frequency  is  finite,  we  know  that  the  random  solution  can 
be  expressed  by  a  Fourier  transform  in  the  time  direction 

M/2  ■  27T 

u(x,t',£)  =  (48) 

n=—M/2 

where  T(£)  is  the  random  period.  It  is  obvious  that  equation  (17)  represents 
the  basic  properties  of  each  random  mode.  We  note  that  gPC  can  effectively 
capture  the  random  behavior  for  some  other  initial  conditions,  e.g.,  Uq{x\  £)  = 
xn.  if  the  polynomial  order  is  large  enough. 

The  long-term  behavior  for  aforementioned  initial  conditions  is  similar  to  the 
spectral  expansion  of  deterministic  functions  with  high  wave  numbers  (see 
[24]).  Since  gPC  is  indeed  a  spectral  expansion  in  terms  of  certain  random 
variables  and  the  time  plays  a  role  similar  to  a  wave  number  in  chaos  expan¬ 
sion,  the  order  of  gPC  must  increase  with  time  to  maintain  a  desired  accuracy 
level. 


3.6  Numerical  results 


Next  we  present  some  numerical  results  for  eM  and  cm-  Let  n  =  1,  v  =  0,  a  =  1 
in  equation  (16).  Due  to  the  periodic  condition  in  physical  space,  we  use  a 
Fourier-collocation  method  to  solve  the  deterministic  PDEs  introduced  by  the 
Galerkin  projection  in  the  gPC  or  ME-gPC  method.  It  has  been  assumed  that 
Qm  <  1  for  the  convergence  of  summation  in  equation  (21).  However,  qm  is  an 
increasing  function  of  t,  which  means  that  the  error  cm  increases  with  time 
and  it  will  reach  0(1)  values  eventually.  In  figure  1,  we  present  the  evolution 
of  the  error  bounds  of  gPC  and  ME-gPC,  respectively.  It  can  be  seen  that  the 
error  of  ME-gPC  increases  at  the  same  speed  as  that  of  gPC.  However,  since 
ME-gPC  is  much  more  accurate  than  gPC,  it  takes  longer  time  for  ME-gPC 
to  reach  0(1)  error.  The  error  cm  of  eighth-order  gPC  is  0(1)  around  t  =  2. 
In  figure  2,  6m  and  the  corresponding  numerical  errors  are  shown,  where  we 
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also  plot  J?+1/2{omrt )  for  comparison.  It  is  known  that  Jf+1  /2(crnnt)  decreases 
exponentially  with  i  when  i  is  much  greater  than  aunt]  otherwise,  there  is  no 
p-convergence.  It  can  be  seen  that  p-convergence  does  not  occur  until  M  >  8 
{qm  <  1)  and  the  rate  of  convergence  is  the  same  as  the  decreasing  rate  of 
J42+i/2(af7r)  when  i  — ■»  oo. 


In  hgure  3  we  demonstrate  theorem  5  numerically.  According  to  theorem  5,  we 
know  that  the  error  of  gPC  at  time  t  should  be  almost  the  same  as  the  error 
of  ME-gPC  of  the  same  polynomial  order  at  time  Nt  for  a  uniform  mesh.  For 
gPC,  we  re-scale  the  time  by  a  factor  N  while  keeping  the  errors  unchanged.  It 
can  be  seen  that  the  re-scaled  error-time  curve  of  gPC  matches  very  well  with 
the  error-time  curve  of  ME-gPC,  and  it  appears  that  the  errors  of  ME-gPC 
are  always  bounded  by  the  shifted  errors  of  gPC.  We  note  that  the  random 
inputs  are  uniform. 


To  examine  if  the  previous  results  extend  to  other  random  distributions,  we 
consider  Beta  and  Gaussian  distributions.  In  hgure  4,  we  plot  the  errors  of 
gPC  and  ME-gPC  versus  time  for  a  Beta  distribution  Beta( 0, 1)  while  the 
time  for  gPC  is  re-scaled  as  before.  It  is  seen  that  the  two  curves  for  gPC  and 
ME-gPC  agree  with  each  other  very  well.  We  also  compare  gPC  and  ME-gPC 
of  polynomial  order  p  =  3  in  hgure  5  when  £  is  of  Gaussian  distribution.  For 
the  ME-gPC  method,  we  hrst  decompose  the  support  of  the  Gaussian  distri¬ 
bution  into  three  random  elements:  ( — oo,  — 6],  [—6,6]  and  [6,  oo).  We  subse¬ 
quently  decompose  the  middle  element  [—6,  6]  while  keeping  the  tail  elements 
unchanged  since  Pr(£  e  (— oo,  — 6]  or  [6,  oo))  =  1.97  x  10“9.  The  inhuence  of 
tail  elements  can  be  observed  only  in  the  early  stage,  which  is  clearly  shown 
in  the  left  plot.  Starting  from  t  ~  0.5,  the  tail  elements  will  not  affect  the 
accuracy  any  more.  Thus,  we  can  drop  the  tail  elements  when  we  re-scale 
the  error  curves  of  gPC  by  the  element  number  N  due  to  their  negligible  er¬ 
ror  contribution.  In  the  plot  on  the  left,  it  can  be  seen  that  the  error  curve 
re-scaled  by  factor  2.5  matches  the  ME-gPC  results  better  than  the  element 
number  IV  =  4.  The  reason  is  that  Gaussian  distribution  is  non-uniform.  If 
the  middle  element  is  decomposed  to  four  equidistant  ones,  we  know  that 
Pr(/[_6  3]  =  1  or  / r3)6]  =  1)  =  2.7  x  10”3,  which  implies  that  the  error  of  ME- 
gPC  is  mainly  controlled  by  elements  [—3,0]  and  [0,3]  (see  equation  (33)). 
Thus  the  scale  factor  should  be  about  2.  In  the  plot  on  the  right,  we  can  see 
that  for  larger  element  numbers,  the  re-scaled  error  curves  of  gPC  agree  well 
with  the  ME-gPC  results,  which  implies  the  relation  presented  in  theorem  5 
will  be  valid  asymptotically  for  the  Gaussian  distribution. 
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Fig.  1.  Evolution  of  error  estimates  of  gPC  and  ME-gPC.  Here,  £  is  uniform  in 
[-1,11- 


Fig.  2.  Convergence  of  gPC  and  J?2+1^2  at  t  =  2.  Here,  £  is  uniform  in  [—1, 1], 

4  Application:  noisy  flow  past  a  stationary  circular  cylinder 

4-1  Random-frequency  inflow  noise 


The  aforemenioned  issues  are  often  encountered  in  numerical  simulations  of 
unsteady  noisy  flows.  We  now  simulate  the  two-dimensional  noisy  flow  past  a 
circular  cylinder  subject  to  the  following  random  boundary  conditions  at  the 
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Fig.  3.  Comparison  of  gPC  and  ME-gPC  of  the  same  polynomial  order  p  =  3. 
Normalized  numerical  errors  of  the  second  moment  are  used.  The  time  of  gPC  is 
multiplied  by  the  number  N  of  random  elements.  Here,  £  is  of  uniform  distribution 
[/[—!,  1],  Left:  ME-gPC  with  N  =  2;  Right:  ME-gPC  with  N  =  3. 


t 


Fig.  4.  Comparison  of  gPC  and  ME-gPC  of  the  same  polynomial  order  p  =  3. 
Normalized  numerical  errors  of  the  second  moment  are  used.  The  time  of  gPC  is 
multiplied  by  the  number  N  of  random  elements.  Here,  £  is  of  Beta  distribution 
Beta( 0, 1)  on  [—1, 1], 

inflow 

u  —  1  +  crY,  v  =  0,  (49) 

where  Y  is  a  uniform  random  variable  of  zero  mean  and  unit  variance  and 
a  is  a  prescribed  constant  indicating  the  degree  of  perturbation.  For  each 
value  of  Y,  there  exists  a  corresponding  Reynolds  number,  which  determines 
a  unique  vortex  shedding  frequency.  In  other  words,  the  shedding  frequency 
in  the  stochastic  simulation  is  random. 

In  figure  6  the  mesh  used  for  the  discretization  in  physical  space  is  shown. 
Neumann  boundary  conditions  (zero  flux)  are  employed  at  the  outflow  and 
periodic  boundary  conditions  in  the  cross-flow  direction.  The  numerical  for¬ 
mulation  of  gPC  for  the  incompressible  Navier-Stokes  equations  was  presented 
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Error 


Fig.  5.  Comparison  of  gPC  and  ME-gPC  of  the  same  polynomial  order  p  =  3. 
Normalized  numerical  errors  of  the  second  moment  are  used.  The  time  of  gPC  is 
re-scaled  by  a  constant.  Here,  £  is  of  Gaussian  distribution  jV(0, 0.2).  Left:  ME-gPC 
with  N  =  4;  Right:  ME-gPC  with  N  =  8, 16. 

in  [25],  where  spectral/ /ip  element  methods  were  employed  to  solve  the  large 
deterministic  PDE  system  produced  by  the  Galerkin  projection  in  gPC.  The 
Reynolds  number  considered  in  this  work  is  Re  =  100. 

We  first  simulate  a  deterministic  case  with  a  =  0  up  to  t  —  1000  to  obtain 
a  fully  developed  flow  and  then  introduce  10%  noise  at  the  inflow.  We  plot 
the  instantaneous  mean  and  variance  of  lift  coefficient  Cl  in  figure  7,  and 
of  drag  coefficient  C-d  in  figure  8.  It  can  be  seen  that  both  the  mean  and 
variance  of  Cl  given  by  ME-gPC  oscillate  periodically  around  a  constant  value 
with  a  decreasing  amplitude  after  a  short  transient  stage.  This  agrees  with  a 
stochastic  model  of  lift  coefficient  developed  in  [18,26],  where  Cl  is  modeled 
by  a  harmonic  signal  with  a  random  frequency.  Based  on  such  a  model,  the 
mean  of  Cl  goes  to  zero  while  the  variance  asymptotes  a  constant  value.  Good 
agreement  between  gPC  and  ME-gPC  is  observed  only  in  the  transient  stage, 
after  which  gPC  begins  to  diverge.  Similar  trends  are  observed  for  the  drag 
coefficient  Cd-  The  study  in  previous  section  shows  that  the  polynomial  order 
of  gPC  must  increase  at  about  a  constant  rate  (see  equation  (38))  to  maintain 
a  certain  accuracy  if  a  random  frequency  is  involved.  In  figure  9,  we  plot  the 
normalized  relative  errors  of  the  variance  of  the  lift  coefficient  using  the  results 
given  by  ME-gPC  with  N  —  20  and  M  —  8  as  a  reference.  It  can  be  seen  that 
the  errors  of  gPC  increase  quickly  to  0(1).  ME-gPC  with  N  =  20  and  p  —  6 
reaches  an  error  of  O(10“2)  at  t  ~  135.  We  note  that  errors  less  than  10-5  are 
not  shown  because  the  output  data  are  truncated  after  the  fifth  digit. 

In  section  3  we  have  shown  that  the  error  of  ME-gPC  at  a  fixed  time  can  be 
estimated  from  that  of  gPC  of  the  same  polynomial  order  but  shifted  by  a 
factor  N.  Here  we  cannot  use  this  result  directly  because  the  decomposition 
of  random  shedding  frequencies  is  not  necessarily  uniform  although  the  noise 
at  the  inflow  is  uniform.  However,  we  can  estimate  the  scaling  factor  from  the 
simulation  results  for  the  errors  of  gPC  and  ME-gPC  of  sixth  order,  which  is 
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Fig.  6.  Schematic  of  the  domain  for  noisy  flow  past  a  circular  cylinder.  The  size 
of  the  domain  is  [—  15D,25D\  x  [— 9D,9D]  and  the  cylinder  is  at  the  origin  with 
diameter  D  =  1.  The  mesh  consists  of  412  triangular  elements. 

about  12.  Using  this  value  we  know  that  the  error  of  ME-gPC  with  N  —  20 
and  M  =  8  at  tU / D  =  150  should  be  roughly  equal  to  the  error  of  eighth-order 
gPC  at  tU / D  =  150/12,  which  is  O(10~3).  Thus,  ME-gPC  with  IV  =  20  and 
M  —  8  can  provide  accurate  results  in  the  range  tU / D  <  150,  corresponding  to 
about  20  shedding  periods  after  the  transient  stage.  In  contrast,  gPC  provides 
accurate  result  up  to  less  than  two  shedding  periods. 

In  figure  10  the  RMS  of  vorticity  is  plotted.  The  global  structure  is  (approxi¬ 
mately)  symmetric  and  the  values  of  RMS  of  vorticity  are  decreasing  gradually 
from  the  front  stagnation  point,  through  the  boundary  layers,  into  the  wake. 
This  suggests  that  the  vorticity  behind  the  cylinder  should  contain  a  harmonic 
signal  M(x,  Y)  cos(2irfv(Y)t)  with  random  frequencies  fv.  The  RMS  of  such 
a  harmonic  response  will  approach  fY  A2(x,  Y)f(Y)/2dY  as  t  — »  oo,  where 
f(Y )  is  the  PDF  of  Y .  Since  the  flux  of  vorticity  decreases  in  the  x  direction 
due  to  viscous  diffusion,  the  value  of  t4(x,  Y)  should  also  decrease  in  the  x 
direction.  This  explains  qualitatively  why  we  only  observe  decreasing  RMS 
values  of  vorticity  in  the  wake  without  the  von  Karmon  vortex  street. 


4-2  Random- amplitude  inflow  noise 


In  this  section  we  consider  another  noisy  boundary  condition  at  the  inflow 

u  =  1  +  cr£  cos  2-nfint,  v  =  0,  (50) 

where  we  add  a  harmonic  signal  with  a  random  amplitude  into  the  inflow.  We 
use  fin  =  0.75 fs,  where  fs  is  the  vortex  shedding  frequency  at  Re  =  100.  Let 
£  be  a  uniform  random  variable  with  zero  mean  and  unit  variance.  We  set  cr 
to  0.1. 

In  figure  11  we  compare  the  mean  and  variance  of  Cl  given  by  gPC  and  ME- 
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Fig.  7.  Evolution  of  mean  (upper)  and  variance  (lower)  of  lift  coefficient.  Y  is  uniform 
in  [—  v/3,  \/3] .  a  =  0.1. 


gPC  of  the  same  order  M  —  8,  where  IV  =  10  for  the  ME-gPC.  We  see  that 
eighth-order  gPC  can  capture  all  the  statistics  up  to  the  second-order  in  the 
range  oitU/D  <  150  in  contrast  to  the  fast  divergence  of  gPC  for  the  random- 
frequency  noise.  Numerical  experiments  show  that  for  the  boundary  condition 
(50)  the  frequency  of  Cl  is  not  sensitive  to  the  boundary  noise,  where  the 
vortex  shedding  frequency  at  Re  =  100  is  dominant.  Thus,  the  error  of  gPC 
increases  much  slower  than  the  first  case.  A  similar  example  is  noisy  flow  past 
an  oscillating  circular  cylinder  [4],  where  the  frequency  is  also  not  sensitive 
to  the  noise  and  thus  gPC  can  do  a  good  job.  Such  observations  imply  that 
the  presence  of  random  frequencies  can  cause  a  significant  degradation  of  the 
performance  of  gPC,  and  thus  employing  ME-gPC  is  necessary  for  convergent 
results. 
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Fig.  8.  Evolution  of  mean  (upper)  and  variance  (lower)  of  drag  coefficient.  Y  is 
uniform  in  [— \/3,  V3].  a  =  0.1. 


0 


tU/D 


Fig.  9.  Evolution  of  relative  errors  of  instantaneous  variance  of  lift  coefficients  given 
by  gPC  and  ME-gPC.  The  results  obtained  from  ME-gPC  with  N  =  20  and  M  =  8 
are  used  as  a  reference. 
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Fig.  10.  Instantaneous  spatial  distribution  of  RMS  of  vorticity  at  tU/D  =  100. 
N  =  20  and  M  =  8. 
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Fig.  11.  Evolution  of  mean  (left)  and  variance  (right)  of  Cl  given  by  gPC  and 
ME-gPC.  Here,  £  is  a  uniform  random  variable  in  [— \/3,  \/3]  with  a  =  0.1. 

5  Summary 


In  this  work  we  studied  the  long-term  behavior  of  gPC  and  ME-gPC  by  fo¬ 
cusing  on  problems  related  to  random  frequencies.  We  first  analyzed  the  one¬ 
dimensional  advection  equation  with  a  uniform  random  transport  velocity,  for 
which  the  error  estimates  of  gPC  and  ME-gPC  were  derived  for  the  Legendre- 
chaos  expansion.  Based  on  the  error  estimates,  we  found  that  ME-gPC  with  a 
uniform  mesh  in  random  space  can  extend  the  valid  integration  time  of  gPC 
by  a  factor  N,  which  is  the  number  of  random  elements.  Subsequently,  we 
extended  this  relation  to  other  random  distributions  and  verified  it  by  numer¬ 
ical  studies  on  Beta  and  Gaussian  distributions.  We  then  simulated  noisy  flow 
past  a  stationary  circular  cylinder,  where  two  different  boundary  conditions  at 
the  inflow  were  considered.  For  the  random-frequency  noise  at  the  inflow,  the 
vortex  shedding  frequency  is  sensitive  to  the  inflow  condition;  thus  gPC  fails 
to  converge  at  early  integration  times  but  ME-gPC  can  capture  the  correct 
random  behavior  effectively.  For  the  random-amplitude  noise  at  the  inflow,  the 
shedding  frequency  is  not  as  sensitive;  thus  gPC  is  valid  for  a  relative  longer 
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time.  This  flow  problem  provides  good  support  for  our  analysis  in  section  3. 
Although  ME-gPC  can  improve  the  performance  of  gPC  for  problems  related 
to  random  frequencies,  it  also  fails  asymptotically.  To  treat  this,  one  choice 
is  to  increase  adaptively  the  number  of  elements  of  ME-gPC  to  maintain  a 
reasonable  accuracy  in  a  desired  range  of  integration  time.  However,  for  high¬ 
dimensional  random  inputs,  the  effectiveness  of  ME-gPC  will  be  weakened 
since  the  number  of  elements  may  increase  fast.  Thus,  the  long-term  behavior 
of  polynomial  chaos  deserves  further  study. 
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